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Abstract 

The inclusion of electromagnetic effects in unquenched QCD can be accomplished 
using ensembles generated in dynamical simulations with pure QCD provided the 
change in the quark determinant induced by a weak electromagnetic field can be effi- 
ciently computed. A stochastic technique for achieving this in the case of dynamical 
domain wall calculations is described. 



1 Introduction 



Previous work in quenched QCD [1, 2] has estabhshed the possibihty of computing quark 
masses including electromagnetic effects, as well as fine structure of hadron multiplets, by 
explicitly coupling the quark fields to a weak electromagnetic U(l) gauge field superimposed 
on pure QCD gauge fields generated in a conventional Monte Carlo simulation. Although the 
pattern of fine structure revealed in these quenched calculations qualitatively matches the 
experimental splittings for the meson and baryon octets, a credible quantitative computation 
clearly requires a fully unquenched treatment of both the gluon and photon dynamics. In this 
paper we present a technique for including electromagnetic effects in preexisting ensembles of 
pure QCD gaugefields generated with dynamical domain wall quarks. Although we describe 
the method for domain wall quarks (preferred in this case for the enhanced chiral symmetry 
and control over quark mass renormalization available in this formulation), the technique 
applies just as well to other fermion discretizations such as Wilson quarks (with or without 
clover improvement). 

To illustrate the basic idea, we note that the path integral for quarks coupled to both 
SU(3) gluons (represented by conventional compact link variables U, with action jS'gi) and an 
abelian photon field A (with a noncompact action 5'em(^) = 4^ J2niJ,u{^ tJ,^nu — Vj/74„^)^,[1], 
in Coulomb gauge) can be written (somewhat schematically): 

Zsu(3).u(i) = / PC/PA det(M(C/,A))e-^«'(^)-^-(^) (1) 

= [VU det(M(^,A = 0))e-^«'(^) [ VA ^ ^^f i^t(^^^'^)) -SMA) ^3) 
J V V , J det(M([/,A = 0)) ^ ' 

where M{U, A) is the matrix defining the lattice quark action in the presence of the SU(3) 
field U together with the abehan photon field A. Evidently, at least in principle, the sim- 
ulation of physical observables in the presence of both QCD and electromagnetism can be 



effected by superimposing photon fields generated with weight e 



Sem(A)+ln 



dct(M(;7,A)) 

det(M(£/,o)) on a pre- 



existing ensemble U generated by a pure SU(3) dynamical QCD simulation. This approach 
will of course only be feasible provided: 



2. there is a reasonably efficient strategy for computing the above determinant ratios for 
a given gluon field U and a suitable ensemble of weak photon fields A. 

Given the extremely expensive nature of full dynamical QCD calculations, especially in 
chirally invariant formulations such as with domain wall or overlap fermions, it is unrealistic 
to expect that dynamical configurations with explicitly simulated photon fields, at a variety 
of quark masses and electric charges, will be available in the foreseeable future. Consequently, 
a procedure which allows us to exploit dynamical configurations generated for electrically 
neutral quarks to study issues of electromagnetic fine structure would clearly be very useful. 

In Section 2, we show how the determinant ratios appearing in (2) may be computed 
exactly, at least for small lattices, for the domain wall formulation of SU(3)xU(l) gauge 
theory, using a Lanczos algorithm. Although this approach is not feasible for large lattices, 
it gives an exact result on smaller lattices which is extremely useful in checking the correctness 
and accuracy of the stochastic approach to computing determinant ratios described in Section 
3. In Section 4 we give further details regarding the computational effort required to extract 
determinant ratios at a given level of accuracy. 

2 Exact Computation of Domain Wall Determinants 

In the domain wall formulation [3, 4] of lattice quark fields, the matrix D{U) defining the 
quadratic form for the fermionic action can be rendered hermitian by premultiplication with 



1 



The variance of the determinant ratios In ( 



det{M{U,A)) 
det(M([/,0)) 



) is tolerably small, and 



the four-dimensional 75 matrix and the time inversion operator R in the fifth dimension 
(assumed of extent L5 henceforth), yielding a matrix M{U) = R'^^D{U) which as an -L5XL5 
dimensional block matrix takes the form (for ease of display, taking ^5=6) 
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where P± = |(1±75), niq is the bare quark mass, 05 the lattice spacing in the fifth dimension, 
H the four-dimensional hermitian Wilson operator with an appropriately chosen domain wall 
mass (or hopping parameter k). It is easy to establish the following trace identities 

Tr(M) = (3) 
IV(M2) = i2V{L,{l-^y + L,{l + 4al)+ml-l) (4) 

where V is the volume of the four-dimensional lattice. These trace identities serve as useful 
checks on the accuracy of the computed spectrum of M: identities involving higher powers 
of M can also be derived, involving plaquette (or even higher loop) sums, but will not be 
needed here. Of course, det(M) = det{D). 

For not too large lattices, the Lanczos algorithm can be used to extract the complete 
spectrum of the hermitian matrix M: for a review of the issues involved in computing 
complete spectra with Lanczos, we refer the reader to [5]. In this paper we shall consider 
5-dimensional lattices dimensioned L^xLs with L =4,6 and L5—I2. For the case L —A, the 
matrix M has dimension 12\^L5=36864, and resolving the complete spectrum to 9 significant 
figures requires on the order of 200,000 Lanczos sweeps. One then finds that Tr(M) ~10~^ 



and the left and right hand sides of (4) agree to 9 significant figures. The computational 
effort required grows as the square of the dimension of M, however, and on larger lattices 
one eventually encounters missed eigenvalues (due to accidental close degeneracies) which 
make this brute-force approach to extracting a full spectrum, and determinant, ultimately 
impractical. However, having an "exact" result for the domain-wall determinant (at least, 
to 9 significant figures) will turn out to be a useful check on the accuracy of the stochastic 
methods described below. 

3 Stochastic Procedure for Determinant Ratios 

The stochastic formalism developed by Golub and coworkers [6] provides an approach to 
the calculation of determinant ratios which allows us to go to larger, and more physically 
realistic lattices. This technique was first applied to the quark determinant problem in lattice 
QCD by Irving and Sexton [7]. In this approach, given a positive definite hermitian matrix 
P, the close connection between Lanczos recursion and Gaussian integration is exploited to 
compute an arbitrary diagonal matrix element < v\f{P)\v > of any differentiable function 
f{P) of the matrix P. Any such matrix element can be written as a spectral sum which 
can then be interpreted as a Riemann-Stieltjes integral computable via the usual variety of 
Gaussian quadrature rules (Gauss, Gauss- Radau, Gauss-Lobatto, etc.). It turns out that 
the appropriate weights for the Gaussian quadrature are directly related to the spectrum 
of the tridiagonal matrix generated by Lanczos recursion on P (for further details, see [6]). 
With a technique for computing diagonal matrix elements in hand, this approach then leads 
immediately to a stochastic method for unbiased estimates of the Tr(/(P)) as the average 
over an ensemble of diagonal expectations of f{P) with respect to a set of random vectors 



Zi > where each vector has components ±1 chosen at random: 



Tr(/(P)) 



E{< z\f{P)\z>) 



(5) 



var{<z\f{P)\z>) = 2El/(^)..r (6) 



The above method can be fruitfully applied to our problem- the evaluation of 




as follows. Choosing 



P = M{U, 0)-^M{U, AfM{U, 0) 



-1 



(7) 



we then have a positive definite hermitian matrix P as required for the algorithm described 
above (for the work described in this paper, we have performed the required inversions of 
M(C/, 0) using the minimum residual algorithm [8]). Choosing f{P) — ln(P), we obtain 
the (log of the square of the) desired determinant ratio. Moreover, for sufficiently weak 
electromagnetic fields A it is clear that P is close to the identity, which has the following 
salutary consequences: 

1. As we shall see shortly, the Lanczos recursion to compute the individual diagonal 
matrix elements < 2;|ln(P)|2; > converges extremely rapidly: in fact 4 or 5 Lanczos 
sweeps (corresponding to 8 or 10 inversions of the domain wall matrix M) are typically 
sufficient to obtain the matrix element to single precision accuracy. 

2. As P is close to the identity, (6) implies that the variance of the diagonal matrix 
elements should also be small, allowing us to obtain a good estimate of the trace with 
an ensemble of manageable size. 

We have studied these issues on a variety of small lattices, generating pure SU(3) 
quenched configurations (at /3=6.0) as the strong interaction background field, and Coulomb 
gauge photon field configurations corresponding to electric charge = —0.1 (roughly, the 



down quark charge). On a 4^^x12 lattice we then typically find In det(M(t/o)) ~ |Trln (P) ~ 
—5. For a specific diagonal element < z\ In {P)\z > with a typical random vector \z >, the 
convergence of the Lanczos estimates is indicated in Table 1, both for 4^x12 and for 6"^xl2 
lattices. In either case, the convergence to seven or eight significant figures occurs after less 
than 10 Lanczos sweeps, with the rapidity of convergence basically insensitive to the lattice 
volume. The computational effort required for the Lanczos recursion is of course dominated 
by the inversions of the domain wall operator M([/, 0). 

Table 1: Convergence of Lanczos estimates for a diagonal element of In (P) 



Lanczos sweep 


< z\\n{P)\z >,L = 4 


<^|ln(P)|2>,L = 6 
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-11.29817352 


-53.44524316 


3 


-11.54470058 


-55.10595906 


4 


-11.54923739 


-55.14478821 


5 


-11.54931955 


-55.14572594 


6 


-11.54932110 


-55.14575078 


7 


-11.54932113 


-55.14575146 


8 


-11.54932113 


-55.14575148 



The computation of the determinant ratio, or equivalently the trace of InP, requires 
repeating the evaluation of the diagonal matrix element in a random vector \z > many 
times, but this part of the procedure is of course trivially parallelizable as the random vectors 
(and resulting estimates of the trace) are completely independent. We show in Figure 1 the 
rate of convergence of the desired determinant ratio E{< z \ ln(P)|2; >)/2 as a function of 
the ensemble size, up to a final sample size of 1000 random vectors, for three statistically 
independent photon field A superimposed on the same background QCD field U. For the 



200 400 600 800 1000 

ensemble size (# of random vectors |z>) 

Figure 1: Ensemble convergence for log(determinant ratio), varying photon field, 4^^x12 

red curve, the exact (log)ratio of determinants, evaluated by computation of the complete 
spectrum as described in the preceding section, gives -6.1522, while the stochastic result 
(again, for 1000 random vectors) is -6.2384±0.1414. For the blue curve, the comparison is 
-5.3791 to -5.4472±0.1346, and for the black curve -5.0609 to -5.0635±0.1388. On a 4^x12 
lattice, each evaluation of < z\f{P)\z > in (5) takes about 20 seconds on a Xeon 2. 8GHz 
PC. Note that the variation of the log determinant ratio is of order unity between different 
configurations. In Figure 2 the same convergence issues are displayed, this time varying both 
the abelian A and nonabelian U background fields. 

We have also repeated the stochastic determination of determinant ratios for a variety 
of 6^x12 lattices (with the same mass and coupling parameters as above), and Figure 3 
shows the convergence from an ensemble of 4000 estimates for three different U(l) fields 
superimposed on a fixed QCD background. The variation of the log(determinant ratio) from 
one U(l) field to another again is of order unity. The total variance (6) turns out to be linear 




-8 - 



200 400 600 800 1000 

ensemble size (# of random vectors |z>) 

Figure 2: Ensemble convergence for log(determinant ratio), varying photon, gluon fields, 4^x12 

in the volume of the lattice: for L=4 we find a variance per site of 0.0063(2), while for L=6 
the variance per site is 0.0066(5). The computed log(determinant ratio) is also of course an 
extensive quantity proportional to the lattice volume, so we conclude that we can achieve a 
fixed relative error in the log(determinant ratio) evaluation with a fixed ensemble size, or a 
fixed absolute error by increasing the number of stochastic estimates linearly with the lattice 
volume. Thus the absolute error on the log(determinant ratio) achieved for the 6^x12 lattices 
with 4000 estimates is about 0.16, comparable to the accuracy achieved with 1000 estimates 
on the smaller lattices (which are about one fifth the volume). As emphasized previously, 
the evaluations of < z\f{P)\z > in (5) for different random vectors \z > are completely 
independent so this average is trivially parallelizable with 100% efficiency. 
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Figure 3: Ensemble convergence for log(determinant ratio), varying photon field, 6^x12 

4 Computational Issues 

We have described two different routes to computing the electromagnetic variation of the 
quark determinant in unquenched QCD with domain wall quarks: the direct, "brute force" 
evaluation of the complete spectrum of the (hermitian) domain wall operator via Lanczos 
recursion, or a stochastic, approximate evaluation using the Golub-Bai-Fahey method. The 
first approach, in which the Lanczos recursion is carried out in double precision, delivers a re- 
sult for the determinant (provided all eigenvalues can be resolved) which is typically accurate 
to a precision 5 or 6 significant places short of double precision, with the loss of accuracy due 
to the need for distinguishing spurious from real eigenvalues [9]. The computational effort 
scales superficially as V"^ for lattice volume V (where the first factor of volume arises from 
the multiplication by M{U, A) at each Lanczos sweep, and the second from the fact that 
evaluation of the complete spectrum requires a number of Lanczos sweeps greater than and 
roughly proportional to V). However, for larger lattices the eigenvalue density rises to the 



point where double precision Lanczos recursions will be unable to resolve many of the eigen- 
values, and the method simply becomes impractical in the absence of hardware-implemented 
extended precision arithmetic. Moreover, the most efficient techniques for diagonalization of 
the tridiagonal Lanczos matrix (e.g. the QL algorithm with implicit shifts [10]) are intrin- 
sically serial by nature, which makes parallelization of this part of the procedure difficult 
(although less efficient parallel algorithms, based on the Sturm sequence property [11], have 
been used successfully for small lattices) . The stochastic approach described in the previous 
section also scales as V"^ with the lattice volume, if we require a fixed absolute error in the 
log(determinant ratio) , but does not require extended precision arithmetic as we go to larger 
lattices. Moreover, parallelization of this approach is trivial, as pointed out above. For the 
lattices studied here (L=4,6,L5=12) the two methods are comparable in computational effort 
(if we require an error in the stochastic approach at the few percent level) , so we conclude 
that the stochastic approach is probably the appropriate technique for larger systems. 
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